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Abstract 

A mesh objective crack band model is implemented in the generalized method of cells (GMC) 
micromechanics model to predict failure of a composite repeating unit cell (RUC). The micromechanics 
calculations are achieved using the MAC/GMC core engine within the ImMAC suite of micromechanics 
codes, developed at the NASA Glenn Research Center. The microscale RUC is linked to a macroscale 
Abaqus/Standard finite element model using the FEAMAC multiscale framework (included in the 
ImMAC suite). The effects of the relationship between the characteristic length of the finite element and 
the size of the microscale RUC on the total energy dissipation of the multiscale model are investigated. A 
simple 2-D composite square subjected to uniaxial tension is used to demonstrate the effects of scaling the 
dimensions of the RUC such that the length of the sides of the RUC are equal to the characteristic length 
of the finite element. These results are compared to simulations where the size of the RUC is fixed, 
independent of the element size. Simulations are carried out for a variety of mesh densities and element 
shapes, including square and triangular. Results indicate that a consistent size and shape must be used to 
yield preserve energy dissipation across the scales. 

1.0 Introduction 

Often, details at the constituent level drive the failure of a composite material. An effective way to 
incorporate these details into a numerical model is through multiscale analysis. However, careful attention 
must be paid to how failure is modeled at the different scales. If the stress-strain response of the material 
enters the post-peak softening regime and the tangent stiffness tensor looses positive definiteness, the 
results become sensitive to changes in the mesh. Even if mesh objective constitutive models are in place 
at the microscales or macroscales, mesh objectivity is not necessarily maintained across the scales, and 
the resulting multiscale model will still exhibit pathological dependence on the size of the discretizations 
in one or more scales due to the lack of a clearly defined characteristic length(s) that bridges the scales 
(Bazant, 2007). 

A multiscale model is presented here which couples the generalized method of cells (GMC) 
micromechanics model (Paley and Aboudi, 1 992), to a macroscopic finite element method (FEM) model. 
The commercial FEM software Abaqus is used as the macroscale platform (Simulia Corp., 2008). 
Whereas, the openly distributed Integrated Multiscale Micromechanics Analysis Code (ImMAC) suite of 
micromechanics codes, developed by the NASA Glenn Research Center (GRC), is utilized to perform the 
microscale calculations (Bednarcyk and Arnold, 2002a; Bednarcyk and Arnold, 2002b). The 
Micromechanics Analysis Code/Generalized Method of Cells (MAC/GMC) core engine, within the 
ImMAC suite, is used to implement and solve the subscale GMC problem, and the FEAMAC multiscale 
framework (Bednarcyk and Arnold, 2006) facilitates the communication between Abaqus and 
MAC/GMC via Abaqus user subroutines. 

A mesh objective crack band model (Bazant and Oh, 1983) was implemented at the subcell level 
within the GMC micromechanics model. With the crack band model, the negative, post peak, tangent 
slope of the stress-strain response of the material is scaled so that the strain energy release rate (SERR) of 
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the subcell equals the critical fracture toughness of the material when the stress in the subcell has 
diminished to zero. The crack band model at the microscale is used to influence the behavior of the 
macroscale, in the multiscale FEAMAC model. Pineda et al. (2012) previously implemented the crack 
band model within the high fidelity generalized method of cells (HFGMC), and demonstrated the 
objectivity of the FIFGMC solution with respect to the subcell meshes used in the repeating unit cell 
(RUC). By extension, it is assumed that mesh objectivity is assured within the GMC framework at the 
microscale, but this does not warrant mesh objectivity at the macroscale. 

The objective of this work is to investigate the effect of the size and shape of the subscale RUC, 
related to the macroscale finite element, on the preservation of energy in the multiscale model. Thus, a 
simple example problem is presented to demonstrate these effects. A square, two-dimensional (2-D), 
plane strain composite section was loaded in uniaxial tension, and failure was allowed to evolve within 
microscale RUCs that were linked to the macroscale with FEAMAC. The overall SERR of the model was 
calculated when the dimensions of the RUC were scaled such that the length of the sides of the RUC are 
equal to the characteristic length of the macroscale finite element, and the SERR of the multiscale model 
was also calculated when the dimensions of the RUC were fixed to a realistic value. The calculations 
were performed for numerous macroscale FEM mesh densities and shapes (including square and 
triangular). 

2.0 Modeling Progressive Failure at the Microscale 

A micromechanical analysis technique, coined the method of cells (MOC), was developed by Aboudi 
(1991); subsequently, Paley and Aboudi (1992) expanded MOC into the generalized method of cells 
(GMC), and later Aboudi et al. (200 1 ) further increased the robustness of this method with the high 
fidelity generalized method of cells (HFGMC). These methods provide semi-closed form solutions for 
determining global anisotropic composite properties in terms of the constituent materials, as well as the 
full three dimensional (3-D) stresses and strains in each of the constituent subcells. 

2.1 The Generalized Method of Cells 

It is assumed that a continuously reinforced composite microstructure can be represented as a 
collection of doubly-periodic RUCs containing a general number of constituents, as shown in Figure 1. 
The RUC is then discretized into Np><N r rectangular subcells, as exhibited in Figure 2. Each of these 
subcells is occupied by one of the constituent materials of the composite. The number of subcells and the 
number of materials are completely general. 



Figure 1. — Representation of the 
doubly-periodic microstructure 
of a composite material. 
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Figure 2. — Discretization of a doubly-periodic RUC. 
Implementation of a Smeared Crack Band Model. 


Displacement and traction continuity are enforced in an average, integral sense at the subcell 
interfaces of a discretized RUC. These continuity conditions are used to formulate a set of semi-analytical 
linear algebraic equations that are solved for the local subcell strains in terms of globally applied strains 
or stresses. Subsequently, local constitutive laws can be utilized to obtain the local subcell stresses, and a 
constitutive law for the effective, homogenized composite can be formulated. GMC and HFGMC were 
later reformulated reducing the total number of u nk nowns required to define the problem by Pindera and 
Bednarcyk (1997) and Aboudi et al. (2012), respectively. 

2.2 Implementation of a Smeared Crack Band Model 

When materials enter the post-peak strain softening regime (i.e., failure), stress is related to strain 
through a negative tangent stiffness. If the tangent stiffness tensor looses positive definiteness, failure 
localizes to the smallest length scale in the continuum problem (Bazant and Cedolin, 1991). In FEM this 
is a single element, and in GMC this is a single subcell. Thus, the post-peak strain softening energy is 
dissipated over the volume of the discretization entity it localizes to. Since a stress-strain relationship 
prescribes the energy density dissipated during the failure process, the total amount of energy dissipated 
in the local volume decreases if the size of that local volume is reduced, and in the limit the structure 
would theoretically require zero energy to fail. This pathological mesh dependence of numerical 
simulations utilizing constitutive laws containing a negative tangent stiffness that relates stresses to 
strains is well documented in the literature (Bazant and Cedolin, 1979; Pietruszczak and Mroz, 1981; 
deBorst, 1987). 

One way to overcome this mesh dependence is the use of non-local constitutive laws, which utilize 
the strain gradients to introduce the influence of surrounding elements into the stress-strain behavior of a 
failing element, preventing localization (Eringen, 1996; Bazant, 1994). A more tractable method involves 
judiciously scaling the post-peak softening slope of the stress-strain curve. Then, the energy density 
dissipated via failure is a function of the characteristic length within the discretized continuum, and the 
strain energy release rate (SERR) (energy required to create a unit area of new crack surface) is preserved 
despite changes in the size of the mesh. This type of constitutive law is referred to as a crack band, or 
smeared crack, model (Bazant and Oh, 1983; de Borst and Nauta, 1985). 

A variation of the crack band model was implemented in HFGMC by Pineda et al. (2012), and was 
also implemented in GMC for this work. Crack band initiation is determined using a maximum principal 
stress criterion. 
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where a\^ y) is the maximum principal stress, and a^ y) is the mode 1 cohesive strength of the material. 

Upon initiation, the crack band is oriented such that the normal to the crack band n\^ y) is aligned with the 
direction of maximum principle stress. Figure 3 displays a crack band within a discretized continuum, and 
a magnified view of a single subcell shows the crack band of width w c oriented perpendicular to n\^ y) . 
The characteristic length of the subcell l[P y) is taken as the length of a line running across the subcell, 
oriented in the direction of n \^ y) . 

The behavior of the crack band is controlled by a traction t versus separation [w] law, shown in 
Figure 4(a). The traction-separation law prescribes that the SERR of the subcell (given by the area under 
the traction separation curve) is equal to the fracture toughness, or critical SERR (G c ), of the material 
when the traction reaches zero. Since no actual discontinuity is present within the GMC subcell, the 
subcell behaves as a continuum and obeys a stress versus strain constitutive law, displayed in Figure 4(b). 
The dissipated strain energy density W F is the area under the stress-strain curve. Additionally, the stress 
state in a subcell is uniform; therefore, the subcell stresses and tractions on the crack band are equivalent. 
To ensure mesh objectivity, the negative, post-peak softening slope of the stress-strain curve E T is 
adjusted so that, for subcell fly which contains the crack band, 

liPr ) w yr) ) = G (Py) (2) 




Figure 3. — (Left) Crack band embedded in a discretized continuum. 
(Right) Magnified subcell displays crack band orientation within 
subcell, as well as characteristic length of subcell. 




(a) 

Figure 4. — (a) Traction separation law governing behavior of crack 
band. (b). Stress-strain behavior of material. 
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where l{( iy ) is the characteristic length of the subcell, defined in Figure 3, W\P‘ ’ is the total strain energy 
density, which is dependent on l\P y) that is dissipated in the subcell when the stress reaches zero, and 

5 

G\c r) is the fracture toughness of the material occupying the subcell (material property). Equation (2) 
ensures objectivity with respect to changes in the subcell mesh by maintaining that the total SERR, upon 
complete failure (i.e., completely diminished stress), in the subcell is fixed and independent of the size of 
the subcell. For complete details on the implementation of the crack band model, please refer to Pineda et 
al. (2012). 

3.0 FEAMAC Multiscale Framework 

A multiscale framework is utilized to simulate the deformational response of fiber-reinforced 
composite structures by modeling the fiber-matrix architecture as an RUC at the microscale using GMC 
and coupling the microscale to the lamina/laminate level (macrsocale) FEM model. A synergistic 
approach is employed which executes concurrent multiscaling in time, but sequential in length (Sullivan 
and Arnold, 2011). The commercial finite element software, Abaqus 6.10-1 (Simulia Corp., 2008), is used 
as the FEM platform, and the MAC/GMC core micromechanics software (Bednarcyk and Arnold, 2002a; 
Bednarcyk and Arnold, 2002b) is used to perform microscale calculations. The scales are linked using the 
FEAMAC software implementation (Bednarcyk and Arnold, 2006), which utilizes various 
Abaqus/Standard user subroutines. Subsequently, a wrapper was developed (FEAMAC/Explicit) which 
achieves the same multiscaling tasks as FEAMAC, except Abaqus/Explicit is used as the FEM platform 
(Pineda et al., 2009). MAC/GMC, FEAMAC, and FEAMAC/Explicit are available in the openly 
distributed ImMAC suite of micromechanics codes, developed by GRC. The number of subcells and 
constituents within an RUC at the microscale is completely general, and as many, or as few, as needed 
may be used to accurately characterize the micro-architecture of the composite. Moreover, any 
constitutive model can be used in the subcells of the RUC to rerpresent the mechanical or thermal 
response of the constituents. 

A schematic displaying a typical multiscale model using FEAMAC is displayed in Figure 5. The 
integration point strains are applied to the RUC and the local subcell fields are determined using GMC 
(this process is referred to as localizaion). If the subcell material behavior is nonlinear, the local stresses 
and strains are used to calculate the local stiffnesses, plastic strains, thermal strains, and/or state variables 
via the local constitutive law. The RUC is then homogenized and the global stiffnesses, plastic strains, 



Figure 5. — Diagram showing coupling of macroscale FEM and microscale GMC models. 
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thermal strains, and/or state variables are computed. The global stresses at the integration point are then 
calculated using these global, homogenized fields, and the effects of any nonlinear subscale phenomena 
are introduced into the macroscale through changes in the integration point stress state. The global 
stresses, material jacobian, and updated state variables are then supplied to the Abaqus user material 
UMAT subroutine. For complete details on the FEAMAC implementation, the reader is referred to 
Bednarcyk and Arnold (2006). 


4.0 Effect of Subscale Size on Global Energy Dissipation 

A mesh objective continuum damage mechanics (CDM) model for post-peak strain softening 
implemented within the GMC micromechanics framework was presented in Section 2.2. Using such a 
model ensures that the solutions obtained are insensitive to changes in the subcell mesh (Pineda et al., 
2012). Nonetheless, this does not guarantee mesh objectivity if this model is linked to a macroscale model 
via a multiscale framework. To ensure mesh objectivity the SERR of all elements must be preserved upon 
changes in the macroscale mesh. 

Typically, an RUC is considered dimensionless as the periodic boundary conditions are meant to 
simulate an infinite array of the RUC. Elowever, the domain of the integration point that the microscale 
model is coupled to has some finite dimensions. As energy is dissipated at the microscale within the 
RUC, the macroscale also exhibits energy dissipation governed by the stress-strain response of the 
homogenized RUC, and this energy is dissipated over the volume of the integration point domain. 

The dissipated energy density at the microscale W RUC can be calculated as the total energy density 
minus the elastic strain energy density. 


W* uc 



}_ RUC RUC 
2 ij ij 


\v RUC 


( 3 ) 


where £ RUC and a RUC are the global, homogenized RUC strain and stress tensor, respectively, and V RUC 

is the volume of the RUC. Thus, the SERR G RUC can be calculated as the total dissipated energy density 
multiplied by the characteristic length of the RUC. 

g ruc = i RUC w RUC (4; 

If the volume of the RUC remains constant for all simulations, then the SERR also remains constant. 

Noting that, in a multiscale model, the integration point stress and strains are equivalent to the global 
stresses and strains of the homogenized RUC, then the dissipated strain energy density is automatically 
preserved across the scales. 


W RP =W RUC (5) 

where W^' P ' is the dissipated energy density of the integration point. The SERR of the integration point 
domain G LP ‘ can be calculated by multiplying Equation (3) by the characteristic length of the integration 
point domain l l ( J ' . 


G IP - =l I c p W* uc 


(6) 


In multiscale simulations that keep the volume of the RUC fixed, the dissipated strain energy density of 
the RUC W RLjC is constant, and thus, the SERR is directly proportional to the characteristic length of the 
integration point domain to which the failure localizes. Consequently as the mesh is refined, the SERR in 
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the localized domains, and thus the entire structure, changes. The results will never converge and are 
pathologically mesh dependent. 

To enact mesh objectivity in the multiscale model, the RUC must preserve the characteristic length of 
the integration point domain. 


, RUC _ ,I.P. 
l C ~ l C 


(?) 


If the characteristic length is preserved, then Equation (6) becomes 

G IP ■ = G ruc (8) 

and the SERR of the RUC is preserved within the integration point domain after homogenization. Since 
the local SERR release rate in the failing subcells is governed by the critical fracture toughness (a fixed 
material property), the global SERR of the RUC is independent of size. By extension, the SERR of the 
integration point domain is also insensitive to changes in size. 

The challenge in achieving mesh objective multiscale analyses lies in correctly defining the global 

characteristic length l RUC , or l [ P ‘ . Locally within the RUC, the characteristic length can be defined for a 
subcell by the orientation of the crack band in the subcell (see Section 2.2). Conversely, in the 
homogenized RUC, there is no single crack band angle to define the characteristic length. This can be 
assuaged by defining the RUC such that the dimensions are exactly the same as the dimensions of the 
integration point domain. Thus, every length that can be defined in the integration point domain is defined 
equivalently in the RUC, and the characteristic length is automatically preserved across the scales, even if 
it is unknown. 

Unfortunately, in the doubly-periodic GMC micromechanics framework, the shapes of the RUCs are 
restricted to rectangles; however, 2-D finite elements may take the shape of any general quadrilateral. The 
following sections investigate the effects of scaling the size of the RUC to correspond to the volume of 
the FEM integration point domain on the total SERR of the multiscale model. To demonstrate mesh 
obj ectivity, the total SERR of the entire model must be preserved as the FEM mesh is changed. 


4.1 Numerical Model 

A simple example problem is used to demonstrate the effects of scaling the dimensions of the 
microscale RUC appropriately on the overall SERR of the multiscale model. Additionally, the effects of 
using incompatible finite elements and RUC shapes are investigated with the same example problem. 
Figure 6 exhibits the multiscale example problem. At the macroscale a 1- by 1-nim square section is 
modeled using a uniform mesh consisting of square, 2-D, plane strain, reduced integration, CPE4R 
elements with dimensions l e x l e . A global x-y-z coordinate system is employed at the macroscale. The 
microscale domain is modeled using a doubly-periodic, 7 subcells x 7 subcells GMC RUC with overall 
dimensions L x L. To maintain the proper aspect ratio, for all RUC sizes, to represent a circular fiber the 
height and length of the RUC must be equal. The RUC consists of 13 fiber subcells (colored blue) and 36 
matrix subcells (colored green), and is meant to represent a square packing architecture. A local x i -x 2 -x 3 
coordinate system is used where X\ is the axial, fiber direction, and Wx; are the transverse directions. The 
Xi/z-directions are coincident, as are the xx/x-directions, and the x 3 /y-directions. 

The elastic properties of the fiber and matrix constituent subcells were chosen so that the overall 
elastic response of the homogenized RUC represented the elastic behavior of an 1M7/8552 lamina. The 
transversely elastic properties of 1M7/8552 are presented in Table 1 (Camanho et al., 2007), where the x-y 
plane in Figure 6 is the plane of isotropy and E z: is the axial (fiber direction) Young’s modulus, E xx is the 
transverse stiffness, v zx is the axial Poisson’s ratio, and G zx is the axial shear modulus. The matrix subcells 
in the GMC RUC initially behave as an elastic, isotropic material, and the fibers subcells are initially 
elastic and transversely isotropic. The initial, elastic properties of the microscale constituents are given in 
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Table 2, where E" is the matrix Young’s modulus, v'" is the matrix Poisson’s ratio, E( x is the axial fiber 

Young’s modulus in the x r direction (see Figure 6 ), E{ 2 is the transverse stiffness of the fiber, v{ 2 and 

1//3 are the axial and transverse Poisson’s ratios, and G{ 2 is the axial shear stiffness. A fiber volume 

fraction V 1 of 59.1 percent was used (Camanho et al., 2003). The transverse and axial shear moduli of the 
fiber were obtained from Goldberg and Gilat (2003). The elastic properties of the matrix and remaining 
u nk nown elastic properties of the fiber were calibrated such that the elastic stiffness of the homogenized 
RUC correlated to the elastic properties of the 1M7/8552 lamina. 


macroscale: FEM mesh 


microscale: 


|< 1mm *| 7x7 GMC subcell mesh 



l e 

Figure 6. — Multiscale simulation utilizing square elements at the macroscale 


TABLE 1— ELASTIC PROPERTIES OF 
IM7/8552 LAMINA (CAMANHO ET AL.. 2007) 


E zz , GPa. 
E x „ GPa. 


G a , GPa 


171.4 

.9.08 

.0.32 

.5.29 


TABLE 2.— ELASTIC PROPERTIES OF IM7 CARBON FIBER AND 
8552 EPOXY MATRIX CONSTITUENTS USED 
IN MICROSCALE GMC RUC 


8552 Matrix Properties 

E m , GPa 


4.97 

aA. 


0.36 

IM7 Fiber Properties 

E ! . . GPa 


286.5 

A,', . GPa 


12.4 

Wz 


0.29 

»a 


0.29 

G { 2 , GPa 


20.0 
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The bottom edge of the FEM domain was pinned in the y-direction, and the left edge was pinned in 
the x-di recti on. A uniform displacement was applied to the right edge of the FEM model to simulate 
transverse loading of a 1- by 1-mrn composite section that stretches infinitely in the global, axial z- 
direction (plane strain conditions). Due to the high degree of nonlinearity associated with failure, the 
implicit, dynamic FEM solver available in Abaqus/Standard was utilized to solve the macroscale FEM 
problem using the keyword * DYNAMIC, with the option APPLICATION = QUASI-STATIC (Simulia 
Corp., 2008). A nominal density of 1.57 kg/m 3 was used in all the finite elements. The total displacement 
(ranging from 0.0175 to 0.05 mm), was applied over a time of 1000 sec. 

Failure is introduced into the matrix subcells of the microscale RUC, and the effects are transferred to 
the macroscale through an overall reduction of the stiffness of the RUC (due to a local reduction in the 
stiffness of the failing matrix subcells). The crack band model presented in Section 2.2 is used to govern 
failure at the microscale. Only the yellow and red elements in Figure 6 were modeled with multiscale 
FEAMAC elements. The green elements were modeled using the Abaqus elastic material model and the 
properties listed in Table 1. Free edge effects and the accuracy of the local fields surrounding the failure 
path front are influenced by the mesh refinement, to some degree, which may affect the failure path in the 
simulation. Therefore, failure was restricted to the yellow and red elements to simplify the problem as 
much as possible and try to ensure that the energy dissipated in the model is a result of the multiscale 
methodology, not other influences. Moreover, the matrix subcells of the RUC located in the red element 
were given a slightly lower cohesive strength than the matrix of the RUCs within the yellow elements. 
This was imposed to imitate an imperfection and facilitate the localization of initial failure into a single 
element. Without this imperfection, all of the FEAMAC elements would fail simultaneously because the 
stress and strain fields are uniform throughout the FEM mesh. After failure initiates in the RUC located in 
the red element, failure progresses into the RUCs located within the adjacent yellow elements until the 
failure reaches the free edge of the FEM mesh. 

The local mode 1 cohesive strength cr[< , and the mode 1 fracture toughness G\q of the matrix 
subcells were calibrated using a macroscale mesh of 3 1 elements x 3 1 elements. The dimensions of the 
microscale RUC were chosen to equal the dimension of the finite element; i.e., L = l e = 0.0323 mm. The 
local parameters were calibrated so that the global ultimate stress and total SERR were close to 60.3 MPa 
and 0.2774 kJ/m 2 , respectively, which are the transverse strength and mode 1, transverse fracture 
toughness for IM7/8552 reported by Camanho et al. (2007). The mode 1 cohesive strength, and mode 1 
fracture toughness of the matrix that were used in the simulations are displayed in Table 3. It should be 
noted, that the critical strain of the matrix subcells of the RUCs located in the yellow elements were given 
a value that was 5 percent higher than that used in the red element. 


TABLE 3.— CRACK BAND FAILURE PARAMETERS 
USED IN 8552 MATRIX SUBCELLS 


Mode I cohesive strength a ( / r) (Weak) 56.7 MPa (54.0 MPa) 

Mode I fracture toughness G\^ r) 1.594 kJ/m 2 


Simulations were executed using various square meshes, and the total, macroscopic SERR was 
calculated upon complete diminishment of the load carrying capability. Two scenarios were investigated. 
The first involved scaling the dimensions of the microscale RUC, so that 

L = l I c P - (9) 

while retaining V J = 59.1 percent. Since reduced integration elements are used, the characteristic length of 
an element l e c , and of the integration point domain in an element l [ J> ' , are equivalent. 
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( 10 ) 


The characteristic length of a 2-D element is calculated as the square root of the area of the element A e 
(Simulia Corp., 2008). 


r c =J7 ( 11 ) 

Thus, for a square, reduced integration element, the characteristic length of the integration point domain 
within an element is equal to the length of one of the sides of the element. 

/£* = /* (12) 

In the second scenario, the dimension of the RUC was kept fixed at L = 5.98 pm (this was calculated 
assuming a fiber diameter of 5 pm, and retaining a 59. 1 percent fiber volume fraction), while the FEM 
mesh was refined. 

To determine an upper bound on the error introduced by having incompatible element and RUC 
shapes, the effects of using 2-D, plane strain, triangular CPE3 elements on the total SERR of the model 
was also investigated. Figure 7 displays a typical macroscale mesh utilizing triangular elements. The 
dimensions of the RUC were still scaled according to Equation (9); however, the definition of the 
characteristic length of the integration point domain is no longer given by Equation (12). Utilizing 
Equation (11), l' c P ' becomes 


/ 


i.p. 

c 



(13) 


for triangular elements. Additionally, the same simulations were conducted with fixed RUC dimension 
L = 5.98 pm. 


macroscale: FEM mesh 

microscale: 


H 1 mm *| 7x7 GMC subcell mesh 



l e 


Figure 7. — Multiscale simulation utilizing triangular elements at the macroscale. 
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4.2 Results 


The total SERR of the model G was calculated using 


G = 




ightedge 


rrightedge 


- dU nghtedge ■ 


j ^ Rrightedge 


rrightedge 


jjrightedge 


(14) 


where RF nghtedge is the sum of the reaction forces at all the nodes on the right edge of the model, 
jrightedge ^ ] cn gth of the right edge (1 mm for all simulations), and if‘z hted z e is the uniform displacement 
applied to the right edge of the model. Figure 8 displays G as a function of the characteristic length l’ c P 
for the various meshes. The results for the square meshes, shown in Figure 6, are indicated with square 
markers; whereas, the results for the triangular meshes (see Figure 7) are indicated with triangular 
markers. Filled markers indicate that the dimensions of microscale RUCs were scaled to satisfy 
Equation (9); whereas empty markers indicated that the RUC dimensions were fixed, L = 5.98 pm. 

The element length that is closest to the realistic size of the RUC ( L = 5.98 pm) is V' c = 5.99 pm 
(167 elements x 167 elements), and is only 0.17 percent larger than the actual RUC length. Thus, the 
scaled solution obtained with l PP ' = 5.99 pm is considered to be the most accurate solution. For this mesh 
G was calculated to be 0.227 kJ/m 2 ; this value will be used as the baseline SERR for comparison with the 
other simulations. The most erroneous G was calculated, from the solution using the coarsest mesh 

Iq P ‘ = 32.3 pm, to be 0.263 kJ/m 2 . Therefore, for the worst solution, G was 16 percent higher than the 
most accurate solution; although, the element length was 438.7 percent larger. As the mesh was refined, G 
approached 0.227 kJ/m 2 . The maximum G, minimum G and percentage difference between the two Gs 
and corresponding l[: p s are presented in Table 4 for all the square and triangular mesh simulations. 
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Figure 8. — Total SERR of multiscale models as a function of 
characteristic length. 
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TABLE 4.— MAXIMUM AND MINIMUM G, AS WELL AS CORRESPONDING , 

AND PERCENTAGE DIFFERENCE BETWEEN GS AND I 1 /' 



Quad. - scaled 

Quad. - unsealed 

Tri. - scaled 

Tri. - unsealed 

rLP. 
l C » 

|im 

G, 

kJ/m 2 

1 IP - 
l C > 

|im 

G 

kJ/m 2 

jLP. 

l c > 

|im 

G, 

kJ/m 2 

1 IP - 
L c » 

|im 

G, 

kJ/m 2 

Max G 

32.3 

0.263 

32.3 

0.702 

33.7 

0.397 

33.7 

0.866 

Min G 

5.99 

0.227 

5.99 

0.230 

5.94 

0.276 

6.15 

0.281 

% Diff. 

438.7 

16.0 

438.7 

204.8 

467.3 

44.0 

448.0 

208.2 


If the dimensions of the RUC were not scaled, the error observed in the SERR was very significant. 
The coarsest mesh exhibited G = 0.702 kJ/m 2 which is an error of 209.6 percent from the baseline. 
Although, when the characteristic length of the integration point domain and the size of the RUC are 

comparable (i.e., l' c P ' = 5.99 pm), G = 0.230 kJ/nr, which is only difference of 1.6 percent from the 
baseline. However, the SERR calculated for the unsealed simulations dropped suddenly when the element 
size approached the size of the RUC. For example, the average G for the next three largest meshes 

( l^ p ' = 6.14, 7.41, and 8.70 pm) is 0.328 kJ/m 2 . The SERR for the coarsest mesh is 204.8 percent larger 
than the finest mesh if the RUC dimensions are not scaled. 

Clearly, scaling the size of the RUC such that it corresponds to the dimensions of the finite element 
minimizes the error. The variation that is observed in simulations with scaling, as the mesh is refined, can 
be attributed to local effects in the macroscopic FEM model such as an increase in the stress 
concentration ahead of the failure path front, as the mesh is refined. 

The simulations incorporating triangular meshes exhibited an overall larger degree of error (when 
compared to the baseline) in the SERR. This is to be expected, since the shape the element and the shape 
of the RUC are incompatible. The characteristic length, calculated using Equation (13), that was closest to 
the assumed actual RUC length, is 5.94 pm. When the RUC dimensions were scaled, G was calculated to 
be 0.276 kJ/m 2 for this characteristic element length, an error of 21.6 percent from the baseline calculated 
using square elements. The largest G calculated (when RUC scaling was enforced), corresponding to the 
coarsest mesh with l PP ' = 33.7 pm, was 0.397 kJ/m 2 . The difference between the minimum G and 
maximum G is 44.0 percent, nearly three times the difference seen in simulations with square elements. 

For the unsealed simulations utilizing a triangular mesh, the largest G calculated was 0.866 kJ/m 2 

corresponding to the coarsest mesh ( l' c P ' = 33.7 pm), and the smallest G calculated was 0.281 kJ/m 2 , 
corresponding to l l c P ' = 6.15 pm. Note that the finest mesh ( l[ JP = 5.94 pm) produced an SERR of 
0.294 kJ/m 2 which is slightly higher than the SERR obtained from the slightly coarser mesh. The 
minimum error was 23.9 percent; however the difference between the maximum G and the minimum G 
was 208.2 percent. The mismatch between the element shape and the RUC shape introduces some error 
when compared to the baseline square solution; thus, the use of triangular elements in an FEAMAC 
simulation shoidd be avoided. However, the error introduced through refining the mesh can be minimized 
by scaling the RUC dimension by the characteristic length of the element. 


5.0 Conclusions 

A smeared crack band model was implemented within the GMC micromechanics framework. The 
crack band model utilizes the characteristic length of the subcells in the micromechanics RUC to preserve 
the SERR of each subcell, ensuring that it is equal to the critical fracture toughness of the material and is 
insensitive to changes in the subcell mesh. This yields mesh objective microscale failure. 

A platform for linking macroscale finite element analysis to microscale GMC computations, the 
FEAMAC framework developed at GRC, is presented. This framework utilizes standard Abaqus and 
MAC/GMC input decks and data cards to setup a multiscale problem. Communication between the 
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MAC/GMC core engine, also developed at GRC, and Abaqus FEM commercial software is facilitated 
through an Abaqus user material UMAT subroutine. 

Use of the crack band model at the microscale ensures the microscale behavior is objective with 
respect to changes in the GMC subcell mesh; however, it does not guarantee that the macroscopic 
solution will be objective with respect to the finite element mesh in a multiscale simulation. A scaling 
technique was developed for preserving the inherent length scale of the macroscale as the analysis 
proceeds into the microscale. The scaling procedure utilizes the characteristic length of the integration 
point domain within an element to scale the dimensions of the RUC such that the volume of the 
integration point domain and the RUC are equivalent. Plane strain, 2-D, uniaxial tension simulations were 
performed using meshes with square, triangular, and rectangular elements with, and without, scaling the 
RUCs for numerous degrees of mesh refinement. The overall SERR, compared to a baseline model, of 
each simulation is used to gauge the accuracy of the models. The technique works well when the elements 
are nearly square in shape, and alleviates a large degree of error observed when the scaling is not 
performed. Solutions using triangular elements exhibited significant error; however, the results were 
fairly insensitive to mesh refinement. Elowever, changing the aspect ratio of rectangular elements 
significantly affected the accuracy of the solutions, and scaling the RUC by the characteristic length of 
the elements did not appropriately adjust the microscale problem to preserve the SERR as the mesh was 
changed. 
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